set.seed(100)

######################################################
###                     Packages                   ###
######################################################
library(foreign); library(lattice); library(memisc); library(fields); library(apsrtable); library(xtable); library(sm); library(plotrix); library(ggplot2); library(nnet); library(reshape2)

######################################################
###                      Data                      ###
######################################################
elites <- read.csv("elites.csv", header=T, sep=',')
names(elites)

######################################################
###               Auxiliary Functions              ###
######################################################
multiplot <- function(..., plotlist=NULL, cols) {
  require(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # Make the panel
  plotCols = cols                          # Number of columns of plots
  plotRows = ceiling(numPlots/plotCols) # Number of rows needed, calculated from # of cols
  
  # Set up the page
  grid.newpage()
  pushViewport(viewport(layout = grid.layout(plotRows, plotCols)))
  vplayout <- function(x, y)
    viewport(layout.pos.row = x, layout.pos.col = y)
  
  # Make each plot, in the correct location
  for (i in 1:numPlots) {
    curRow = ceiling(i/plotCols)
    curCol = (i-1) %% plotCols + 1
    print(plots[[i]], vp = vplayout(curRow, curCol ))
  }
  
}

######################################################
###               Missing Solutions                ###
######################################################
elites$p71[elites$p71=="NS" | elites$p71=="NR"] <- NA
elites$p72[elites$p72=="NS" | elites$p72=="NR"] <- NA
elites$p32[elites$p32=="NS" | elites$p32=="NR" | elites$p32=="outro"] <- NA
elites$p35[elites$p35==88 | elites$p35==99 | elites$p35==77] <- NA
elites$p57[elites$p57=="NS" | elites$p57=="NR"] <- NA
elites$p70[elites$p70=="NS" | elites$p70=="NR"] <- NA
elites$p69[elites$p69==99] <- NA

### New Data Set ###
state <- as.numeric(elites$uf) #5.BA,6.CE,9.Go,11.MT,13.MG,14.PA,17.PE,19.RJ,21.RS,24.SC,25.SP,27.TO
state.factor <- factor(state, levels= c(5,6,9,11,13,14,17,19,21,24,25,27),
labels = c("BA","CE","GO","MT","MG","PA","PE","RJ","RS","SC","SP","TO")) 
party <- as.numeric(elites$part_recod)
strategy1 <- as.numeric(elites$p71) #1.all leaders, 2.gov leaders, 3.deputy
strategy1[strategy1==4] <- 0 #outra
strategy1[strategy1==1 | strategy1==2] <- 1 #leader
strategy1[strategy1==3] <- 2 #deputy
strategy2 <- as.numeric(elites$p72) #1.pork deputy, 2.cabinet, 3.portfolios deputies
strategy2[strategy2==4] <- 0 #outra
strategy2[strategy2==1] <- 1 #pork
strategy2[strategy2==2 | strategy2==3] <- 2 #portfolios
strategy <- numeric()
strategy[(strategy1==0 & strategy2==0) | (strategy1==1 & strategy2==0) | (strategy1==2 & strategy2==0) | (strategy1==0 & strategy2==1) | (strategy1==0 & strategy2==2)] <- 0 #outra
strategy[strategy1==1 & strategy2==2] <- 4 #portf/leader
strategy[strategy1==2 & strategy2==1] <- 3 #pork/deputy
strategy[strategy1==1 & strategy2==1] <- 2 #pork/leader
strategy[strategy1==2 & strategy2==2] <- 1 #portf/deputy
constituency <- as.numeric(elites$p32) #1.gov,2.opo,3.party,4.region,5.group,6.popul
constituency1 <- numeric()
constituency1[constituency==5] <- 0 # Sector/Groups
constituency1[constituency==4] <- 1 # concentrated
constituency1[constituency==1 | constituency==2 | constituency==3 | constituency==6] <- 2 # others
leader.power <- as.numeric(elites$p57) #1.weak, 2.strong, 3.middle
leader.power2 <- numeric()
leader.power2[leader.power==1] <- 0
leader.power2[leader.power==3] <- 1
leader.power2[leader.power==2] <- 2
ideology.dep <- as.numeric(elites$p34) #1.left and 10.right
gov.opo <- as.numeric(elites$p69) #1.gov, 5.independ and 10.opo
relation.gov.opo <- as.numeric(elites$p70) #1.competi,2.cooperat,3.depends
relation.gov.opo[relation.gov.opo==2] <- 0
relation.gov.opo[relation.gov.opo==1] <- 2
relation.gov.opo[relation.gov.opo==3] <- 1
first.term <- as.numeric(elites$p79) #1.yes,2.no
first.term[first.term==2] <- 0
first.term[first.term==1] <- 1
gender <- as.numeric(elites$p93) #1.male, 2.female
school <- as.numeric(elites$p106) #0.iliteracy and 8.master/phd
constituency.factor <- as.factor(constituency1)
leader.power.fac <- as.factor(leader.power2)
elites2 <- as.data.frame(cbind(strategy,constituency.factor,leader.power2,ideology.dep, gov.opo, relation.gov.opo, first.term, gender, school, state.factor))

elites2$state.factor <- factor(elites2$state.factor)
elites2$constituency.factor <- factor(elites2$constituency.factor)

head(elites2)

######################################################
###                  Descriptive                   ###
######################################################
# Sample #
table(elites$uf)

# Strategy #
table(strategy)
prop.table(table(strategy))
prop.table(table(strategy1))
prop.table(table(strategy2))
counts <- table(strategy, state)
counts <- rbind(counts, apply(X=counts, MARGIN=2, FUN=sum))
m <- matrix(NA, nrow=6, ncol=12)
for(i in 1:12){
  m[,i] <- counts[,i]/counts[6,i]
}
counts <- m[-6,]
colnames(counts) <- c("BA","CE","GO","MT","MG","PA","PE","RJ","RS","SC","SP","TO")
rownames(counts) <- c("Other", "Portf/Deputy", "Pork/Party", "Pork/Deputy", "Portf/Party")

par(mar = c(5, 3, 2, 2))
barplot(counts, xlim=c(0,21), names.arg=c("BA","CE","GO","MT","MG","PA","PE","RJ","RS","SC","SP","TO"), xlab="States", col=c(1,"gray26","gray50","gray83","gray100"), las=2, legend.text=T)

# Constituency #
table(constituency1)
counts2 <- table(constituency1, state)
counts2 <- rbind(counts2, apply(X=counts2, MARGIN=2, FUN=sum))
m <- matrix(NA, nrow=4, ncol=12)
for(i in 1:12){
  m[,i] <- counts2[,i]/counts2[4,i]
}
counts2 <- m[-4,]
colnames(counts2) <- c("BA","CE","GO","MT","MG","PA","PE","RJ","RS","SC","SP","TO")
rownames(counts2) <- c("Groups","Geographical","Universal")

barplot(counts2, xlim=c(0,21), names.arg=c("BA","CE","GO","MT","MG","PA","PE","RJ","RS","SC","SP","TO"), xlab="States", col=c(225,200,0), las=2, legend.text=T)


# Centralization #
table(leader.power2)
counts3 <- table(leader.power2, state)
counts3 <- rbind(counts3, apply(X=counts3, MARGIN=2, FUN=sum))
m <- matrix(NA, nrow=4, ncol=12)
for(i in 1:12){
  m[,i] <- counts3[,i]/counts3[4,i]
}
counts3 <- m[-4,]
colnames(counts3) <- c("BA","CE","GO","MT","MG","PA","PE","RJ","RS","SC","SP","TO")
rownames(counts3) <- c("Weak","Medium","Strong")

barplot(counts3, xlim=c(0,21), names.arg=c("BA","CE","GO","MT","MG","PA","PE","RJ","RS","SC","SP","TO"), xlab="States", col=c(225,200,0), las=2, legend.text=T)


# Ideology #
table(ideology.dep)
ideology.dep2 <- na.omit(ideology.dep)
ideology.dep.to <- na.omit(ideology.dep2[state==27])

par(mfrow=c(4,3))
plot(density(ideology.dep2[state==5]), xlab="Ideology",main="Bahia", xlim=c(0,10), ylim=c(0,0.40))
plot(density(ideology.dep2[state==6]), xlab="Ideology",main="Ceara", xlim=c(0,10), ylim=c(0,0.40))
plot(density(ideology.dep2[state==9]), xlab="Ideology",main="Goias", xlim=c(0,10), ylim=c(0,0.40))
plot(density(ideology.dep2[state==11]), xlab="Ideology",main="Mato Grosso", xlim=c(0,10), ylim=c(0,0.40))
plot(density(ideology.dep2[state==13]), xlab="Ideology",main="Minas Gerais", xlim=c(0,10), ylim=c(0,0.40)) 
plot(density(ideology.dep2[state==14]), xlab="Ideology",main="Para", xlim=c(0,10), ylim=c(0,0.40))
plot(density(ideology.dep2[state==17]), xlab="Ideology",main="Pernambuco", xlim=c(0,10), ylim=c(0,0.40))
plot(density(ideology.dep2[state==19]), xlab="Ideology",main="Rio de Janeiro", xlim=c(0,10), ylim=c(0,0.40)) 
plot(density(ideology.dep2[state==21]), xlab="Ideology",main="Rio Grande do Sul", xlim=c(0,10), ylim=c(0,0.40))
plot(density(ideology.dep2[state==24]), xlab="Ideology",main="Santa Catarina", xlim=c(0,10), ylim=c(0,0.40))
plot(density(ideology.dep2[state==25]), xlab="Ideology",main="Sao Paulo", xlim=c(0,10), ylim=c(0,0.40))
plot(density(ideology.dep.to), xlab="Ideology",main="Tocantins", xlim=c(0,10), ylim=c(0,0.40))

par(mfrow=c(1,1))

#table2 <- aggregate(mean(ideology.dep, na.omit=T)~party)
#colnames(table2) <- c("Party","Ideology")
#table2 <- sort(table2,by=~Ideology)
#rownames(table2) <- c("PT","PPE","PDT","PSB","PPS","PSDB","PPC","SP","PMDB","PSC","PV","PTB","PP","PPD","DEM","PR")
#plot(table2[c(1,3,4,5,6,9,10,11,12,13,15,16),2],rep(1,12),axes=F,xlab="Ideology",ylab="")
#abline(h=1)
#text(table2[c(1,3,4,5,6,9,10,11,12,13,15,16),2],rep(1,12), c("PT","PDT","PSB","PPS","PSDB","PMDB","PSC","PV","PTB","PP","DEM","PR"), cex=0.6, pos=c(1,1,3,1,1,3,1,3,1,1,3,1), col="red")

summary.stats <- xtable(t(stats(elites2)))
summary.stats

######################################################
###                    Run Model                   ###
######################################################
res1 <- multinom(strategy~factor(constituency.factor),data=elites2, Hess=T)
summary(res1)

res2 <- multinom(strategy~factor(constituency.factor)+factor(leader.power.fac),data=elites2,Hess=T)
summary(res2)

res3 <- multinom(strategy~factor(constituency.factor)+factor(leader.power.fac)+ideology.dep,data=elites2,Hess=T)
summary(res3)

res4 <- multinom(strategy ~ state.factor + constituency.factor + leader.power2 + ideology.dep + gov.opo + relation.gov.opo + first.term + gender + school - 1, data=elites2, Hess=T)
summary(res4)


######################################################
###                    Analysis                    ###
######################################################
### Non-parametric bootstrap for strategies pork.deputy | portf/party ###
simus <- 100 # number of simulations
N <- dim(elites2)[1]
npb.par <- matrix(NA,simus,84)

for (t in 1:simus) {
  thisdta <- elites2[sample(1:N,N,replace=TRUE),]
  res <- multinom(strategy ~ state.factor + constituency.factor + leader.power2 + ideology.dep + gov.opo + relation.gov.opo + first.term + gender + school - 1, data=thisdta, Hess=T)
  npb.par[t,] <- as.vector(t(summary(res)$coefficients))
}

# In npb.par row is each regression, and column each parameter estimated
# Results of Non-parametric Bootstrap
processres <- function(simres) {
  z <- t(apply(simres, 2, function(x) c(mean(x), median(x), sd(x), quantile(x, c(0.05, 0.95)))))
  rownames(z) <- colnames(simres)
  colnames(z) <- c("mean", "median", "sd", "5%", "95%")
  z
}

# Results
processres(npb.par)

# Pork.dep/portf.party #
beta1.strat <- npb.par[,55] - npb.par[,76] # concentrated cons
mean(beta1.strat); sd(beta1.strat); quantile(beta1.strat,c(0.05,0.95))

beta2.strat <- npb.par[,56] - npb.par[,77] # scattered cons
mean(beta2.strat); sd(beta2.strat); quantile(beta2.strat,c(0.05,0.95))

beta3.strat <- npb.par[,57] - npb.par[,78] # party centra
mean(beta3.strat); sd(beta3.strat); quantile(beta3.strat,c(0.05,0.95))

beta4.strat <- npb.par[,58] - npb.par[,79] # ideology
mean(beta4.strat); sd(beta4.strat); quantile(beta4.strat,c(0.05,0.95))

################################################################
###           Marginalized Predicted Probabilities           ###  
################################################################
###      Concentrated Constituency + Decentralized Party     ###
################################################################
nd <- data.frame(state.factor=factor('5', levels=levels(elites2$state.factor)), constituency.factor=factor('2', levels=levels(elites2$constituency.factor)), leader.power2=0, ideology.dep=c(1,2,3,4,5,6,7,8,9,10), gov.opo=mean(gov.opo, na.rm=TRUE), relation.gov.opo=mean(relation.gov.opo, na.rm=T), first.term=mean(first.term), gender=mean(gender), school=mean(school))

## replicate newdataset setting state equal to 1-12 and store 12 datasets in a list
newdat <- rep(list(nd), 12)
newdat <- lapply(seq_along(newdat), function(i) {x <- newdat[[i]]; x$state.factor[] <- i; return(x)})

## get predicted probabilities for every state level
p1 <- lapply(newdat, function(x) predict(res4, type = "probs", newdata=x))

## weight each dataset by the proportion of states
(weights <- with(elites, table(state.factor)/length(state.factor)))
p1 <- lapply(seq_along(p1), function(i) p1[[i]] * weights[i])

## calculate the mean, I can use addition because I used the proportions, which sum to 1
finalp1 <- Reduce(`+`, p1)
test <- cbind(nd, finalp1)
longtest <- melt(test, measure.vars = as.character(0:4), value.name="Probability", variable.name = "outcome")

plot1 <- ggplot(longtest, aes(x = ideology.dep, y = Probability, color = outcome)) + geom_line(aes(linetype=outcome), size=3) + labs(title = "Marginalized Predicted Probabilities \n Concentrated Constituency + Decentralized Party") + ylab("Probability") + xlab("Ideology") + ylim(0,.55) + theme_bw() + theme(axis.title.x = element_text(size = 20, vjust = .25)) + theme(axis.title.y = element_text(size = 20, angle=90)) + theme(legend.position="bottom") + theme(strip.text.x = element_text(size=20, face="bold")) + scale_color_discrete(name='', breaks=c('0', '1', '2', '3', '4'), labels=c("Other", "Port/Dep", "Pork/Leader", 'Pork/Dep', 'Port/Leader')) + theme(legend.text = element_text(size = 16, face = "bold"))

################################################################ 
###        Scattered Constituency + Centralized Party        ###
################################################################
nd <- data.frame(state.factor=factor('1', levels=levels(elites2$state.factor)), constituency.factor=factor('3', levels=levels(elites2$constituency.factor)), leader.power2=2, ideology.dep=c(1,2,3,4,5,6,7,8,9,10), gov.opo=mean(gov.opo, na.rm=TRUE), relation.gov.opo=mean(relation.gov.opo, na.rm=T), first.term=mean(first.term), gender=mean(gender), school=mean(school))

## replicate newdataset setting state equal to 1-12 and store 12 datasets in a list
newdat <- rep(list(nd), 12)
newdat <- lapply(seq_along(newdat), function(i) {x <- newdat[[i]]; x$state.factor[] <- i; return(x)})

## get predicted probabilities for every state level
p1 <- lapply(newdat, function(x) predict(res4, type = "probs", newdata=x))

## weight each dataset by the proportion of states
(weights <- with(elites, table(state.factor)/length(state.factor)))
p1 <- lapply(seq_along(p1), function(i) p1[[i]] * weights[i])

## calculate the mean, I can use addition because I used the proportions, which sum to 1
finalp1 <- Reduce(`+`, p1)
test <- cbind(nd, finalp1)
longtest <- melt(test, measure.vars = as.character(0:4), value.name="Probability", variable.name = "outcome")

plot2 <- ggplot(longtest, aes(x = ideology.dep, y = Probability, colour = outcome)) + geom_line(aes(linetype=outcome),size=3) + labs(title = "Marginalized Predicted Probabilities \n Scattered Constituency + Centralized Party") + ylab("Probability") + xlab("Ideology") + ylim(0,.55) + theme_bw() + theme(axis.title.x = element_text(size = 20, vjust = .25)) + theme(axis.title.y = element_text(size = 20, angle=90)) + theme(legend.position="bottom") + theme(strip.text.x = element_text(size=20, face="bold")) + scale_color_discrete(name='', breaks=c('0', '1', '2', '3', '4'), labels=c("Other", "Port/Dep", "Pork/Leader", 'Pork/Dep', 'Port/Leader')) + theme(legend.text = element_text(size = 16, face = "bold"))

multiplot(plot1, plot2, cols=2)

################################################################ 
###        Scattered Constituency + Decentralized Party      ###
################################################################
nd <- data.frame(state.factor=factor('1', levels=levels(elites2$state.factor)), constituency.factor=factor('3', levels=levels(elites2$constituency.factor)), leader.power2=0, ideology.dep=c(1,2,3,4,5,6,7,8,9,10), gov.opo=mean(gov.opo, na.rm=TRUE), relation.gov.opo=mean(relation.gov.opo, na.rm=T), first.term=mean(first.term), gender=mean(gender), school=mean(school))

## replicate newdataset setting state equal to 1-12 and store 12 datasets in a list
newdat <- rep(list(nd), 12)
newdat <- lapply(seq_along(newdat), function(i) {x <- newdat[[i]]; x$state.factor[] <- i; return(x)})

## get predicted probabilities for every state level
p1 <- lapply(newdat, function(x) predict(res4, type = "probs", newdata=x))

## weight each dataset by the proportion of states
(weights <- with(elites, table(state.factor)/length(state.factor)))
p1 <- lapply(seq_along(p1), function(i) p1[[i]] * weights[i])

## calculate the mean, I can use addition because I used the proportions, which sum to 1
finalp1 <- Reduce(`+`, p1)
test <- cbind(nd, finalp1)
longtest <- melt(test, measure.vars = as.character(0:4), value.name="Probability", variable.name = "outcome")

plot3 <- ggplot(longtest, aes(x = ideology.dep, y = Probability, colour = outcome)) + geom_line(aes(linetype=outcome), size=3) + labs(title = "Marginalized Predicted Probabilities \n Scattered Constituency + Decentralized Party") + ylab("Probability") + xlab("Ideology") + ylim(0,.55) + theme_bw() + theme(axis.title.x = element_text(size = 18, vjust = .25)) + theme(axis.title.y = element_text(size = 18, angle=90)) + theme(legend.position="bottom") + theme(strip.text.x = element_text(size=20, face="bold")) + scale_color_discrete(name='', breaks=c('0', '1', '2', '3', '4'), labels=c("Other", "Port/Dep", "Pork/Leader", 'Pork/Dep', 'Port/Leader')) + theme(legend.text = element_text(size = 16, face = "bold"))


################################################################ 
###      Concentrated Constituency + Centralized Party       ###
################################################################
nd <- data.frame(state.factor=factor('1', levels=levels(elites2$state.factor)), constituency.factor=factor('2', levels=levels(elites2$constituency.factor)), leader.power2=2, ideology.dep=c(1,2,3,4,5,6,7,8,9,10), gov.opo=mean(gov.opo, na.rm=TRUE), relation.gov.opo=mean(relation.gov.opo, na.rm=T), first.term=mean(first.term), gender=mean(gender), school=mean(school))

## replicate newdataset setting state equal to 1-12 and store 12 datasets in a list
newdat <- rep(list(nd), 12)
newdat <- lapply(seq_along(newdat), function(i) {x <- newdat[[i]]; x$state.factor[] <- i; return(x)})

## get predicted probabilities for every state level
p1 <- lapply(newdat, function(x) predict(res4, type = "probs", newdata=x))

## weight each dataset by the proportion of states
(weights <- with(elites, table(state.factor)/length(state.factor)))
p1 <- lapply(seq_along(p1), function(i) p1[[i]] * weights[i])

## calculate the mean, I can use addition because I used the proportions, which sum to 1
finalp1 <- Reduce(`+`, p1)
test <- cbind(nd, finalp1)
longtest <- melt(test, measure.vars = as.character(0:4), value.name="Probability", variable.name = "outcome")

plot4 <- ggplot(longtest, aes(x = ideology.dep, y = Probability, colour = outcome)) + geom_line(aes(linetype=outcome), size=3) + labs(title = "Marginalized Predicted Probabilities \n Concentrated Constituency + Centralized Party") + ylab("Probability") + xlab("Ideology") + ylim(0,.55) + theme_bw() + theme(axis.title.x = element_text(size = 18, vjust = .25)) + theme(axis.title.y = element_text(size = 18, angle=90)) + theme(legend.position="bottom") + theme(strip.text.x = element_text(size=20, face="bold")) + scale_color_discrete(name='', breaks=c('0', '1', '2', '3', '4'), labels=c("Other", "Port/Dep", "Pork/Leader", 'Pork/Dep', 'Port/Leader')) + theme(legend.text = element_text(size = 16, face = "bold"))

multiplot(plot1, plot2, plot3, plot4, cols=2)


##############################################################
###                  Predicted Probability                 ### 
##############################################################
# Concentrated Constituency + Centralized Party + left
nd <- data.frame(state.factor=factor('1', levels=levels(elites2$state.factor)), constituency.factor=factor('2', levels=levels(elites2$constituency.factor)), leader.power2=2, ideology.dep=c(1), gov.opo=mean(gov.opo, na.rm=TRUE), relation.gov.opo=mean(relation.gov.opo, na.rm=T), first.term=mean(first.term), gender=mean(gender), school=mean(school))

## replicate newdataset setting state equal to 1-12 and store 12 datasets in a list
newdat <- rep(list(nd), 12)
newdat <- lapply(seq_along(newdat), function(i) {x <- newdat[[i]]; x$state.factor[] <- i; return(x)})

## get predicted probabilities for every state level
p1 <- lapply(newdat, function(x) predict(res4, type = "probs", newdata=x))

## weight each dataset by the proportion of states
(weights <- with(elites, table(state.factor)/length(state.factor)))
p1 <- lapply(seq_along(p1), function(i) p1[[i]] * weights[i])

## calculate the mean, I can use addition because I used the proportions, which sum to 1
finalp1 <- Reduce(`+`, p1)

# Concentrated Constituency + Decentralized Party + left
nd <- data.frame(state.factor=factor('1', levels=levels(elites2$state.factor)), constituency.factor=factor('2', levels=levels(elites2$constituency.factor)), leader.power2=0, ideology.dep=c(1), gov.opo=mean(gov.opo, na.rm=TRUE), relation.gov.opo=mean(relation.gov.opo, na.rm=T), first.term=mean(first.term), gender=mean(gender), school=mean(school))

## replicate newdataset setting state equal to 1-12 and store 12 datasets in a list
newdat <- rep(list(nd), 12)
newdat <- lapply(seq_along(newdat), function(i) {x <- newdat[[i]]; x$state.factor[] <- i; return(x)})

## get predicted probabilities for every state level
p1 <- lapply(newdat, function(x) predict(res4, type = "probs", newdata=x))

## weight each dataset by the proportion of states
(weights <- with(elites, table(state.factor)/length(state.factor)))
p1 <- lapply(seq_along(p1), function(i) p1[[i]] * weights[i])

## calculate the mean, I can use addition because I used the proportions, which sum to 1
finalp2 <- Reduce(`+`, p1)

# Scattered Constituency + Centralized Party + left
nd <- data.frame(state.factor=factor('1', levels=levels(elites2$state.factor)), constituency.factor=factor('3', levels=levels(elites2$constituency.factor)), leader.power2=2, ideology.dep=c(1), gov.opo=mean(gov.opo, na.rm=TRUE), relation.gov.opo=mean(relation.gov.opo, na.rm=T), first.term=mean(first.term), gender=mean(gender), school=mean(school))

## replicate newdataset setting state equal to 1-12 and store 12 datasets in a list
newdat <- rep(list(nd), 12)
newdat <- lapply(seq_along(newdat), function(i) {x <- newdat[[i]]; x$state.factor[] <- i; return(x)})

## get predicted probabilities for every state level
p1 <- lapply(newdat, function(x) predict(res4, type = "probs", newdata=x))

## weight each dataset by the proportion of states
(weights <- with(elites, table(state.factor)/length(state.factor)))
p1 <- lapply(seq_along(p1), function(i) p1[[i]] * weights[i])

## calculate the mean, I can use addition because I used the proportions, which sum to 1
finalp3 <- Reduce(`+`, p1)

# Scattered Constituency + Decentralized Party + left
nd <- data.frame(state.factor=factor('1', levels=levels(elites2$state.factor)), constituency.factor=factor('3', levels=levels(elites2$constituency.factor)), leader.power2=0, ideology.dep=c(1), gov.opo=mean(gov.opo, na.rm=TRUE), relation.gov.opo=mean(relation.gov.opo, na.rm=T), first.term=mean(first.term), gender=mean(gender), school=mean(school))

## replicate newdataset setting state equal to 1-12 and store 12 datasets in a list
newdat <- rep(list(nd), 12)
newdat <- lapply(seq_along(newdat), function(i) {x <- newdat[[i]]; x$state.factor[] <- i; return(x)})

## get predicted probabilities for every state level
p1 <- lapply(newdat, function(x) predict(res4, type = "probs", newdata=x))

## weight each dataset by the proportion of states
(weights <- with(elites, table(state.factor)/length(state.factor)))
p1 <- lapply(seq_along(p1), function(i) p1[[i]] * weights[i])

## calculate the mean, I can use addition because I used the proportions, which sum to 1
finalp4 <- Reduce(`+`, p1)

## Table that aggregate result
tb <- cbind(c('Concent/Central', 'Concent/Decent', 'Scatter/Central', 'Scatter/Decent'), rbind(round(finalp1, digits=2), round(finalp2, digits=2), round(finalp3, digits=2), round(finalp4, digits=2)))
t(tb)

tb2 <- cbind(c('Concent/Decent', 'Scatter/Central'), rbind(round(finalp2, digits=2), round(finalp3, digits=2)))
colnames(tb2) <- c('', "Other", "Port/Dep", "Pork/Leader", 'Pork/Dep', 'Port/Leader')
t(tb2)